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ABSTRACT 

In paper I of this series, we examined triaxial collapse in terms of the dynamics of 
eigenvalues of three important tensors: the Hessian of the gravitational potential, the 
tensor of velocity derivatives and the deformation tensor. The first paper focussed 
on the joint gravity-velocity dynamics and here we focus on the deformation tensor, 
which is directly related to the axes’ evolution. We examine the evolution of the 
minor to major and intermediate to major axes ratios (s and q) and the triaxiality 
parameter T as function of mass scale and redshift. We find that the ellipticity and 
prolateness increase with decreasing mass scale and decreasing redshift. These trends, 
while in agreement with previous analytic studies, contradict numerical simulations. 
Nevertheless, we find that a suitable transformation of s, motivated by the scaling 
used in recent analysis of the Millennium XXL simulations by Bonamigo et al (2014), 
has a universal log-normal distribution function that matches their numerical results. 
Similarly, the transformation q = (q — s)/( 1 — s) also has a universal beta distribution 
that is valid over a decade in mass range and over a wide range of redshift scales, 
indicating that the variable q can be thought of as an invariant of the phase space 
dynamics. 
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1 INTRODUCTION 


Observations have established beyond any doubt that dark matter haloes are well described by a triaxial geometry. The 
evidence has come from a variety of probes including the optical and X-ray surface brightness, indirect observations though 


the Sunyaev-Zeldovich effect and weak and strong gravitational lensing (see review by 

Limousin et al. 2013 and references 

therein). Numerical simulations have also confirmed triaxiality (e.g., Jing & Suto. 2002 

Allgood et al. 2006 Schneider et al. 

2012) and ellipsoidal halo-shape finders have proven to give more realistic haloes than their spherical counterparts (Despali, 


currently accepted hierarchical model of structure formation small scale structures collapse first and larger structures are 
formed through mass accretion and major mergers. The details of the halo shape distribution provide useful information 
regarding this process of structure formation and evolution. The second major reason is from the point of view of precision 
cosmology. For example, understanding the systematics due to intrinsic ellipticity correlations is a crucial part of any weak 


lensing measurement (e.g., Refregier 2003 Oguri & Keeton 


2004 


Joachimi et al. 20131 and important in determining the 


Hubble constant from cluster shapes (Wang fc Fan|2006 Kawahara et al.|2008|. Cluster and void ellipticities have also been 


proposed as a probe to constrain cosmological parameters such as os and the dark energy equation of state (Hoi et al. 2006; 
Lee|2006 Park fe Lee|2007 Lavaux fe Wandelt||2010 1. 


Over the years, numerical simulations have become increasingly sophisticated and efficient at the task of modeling dark 


matter halo shapes (Dubinski & Carlberg||1991 Bullock||2002 Jing Sz Suto 2002 Kasun &; Evrard 

2005 Hopkins, Bahcall, 

& Bode 2005 Allgood et al.||2006 Schneider, Frenk, & Cole||2012 Despali, Giocoli, & Tormen||2014 

Bonamigo et al. 2014). 
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2 Sharvari Nadkarni-Ghosh and Akshat Singhal 


Using this information for precision cosmology involves exploring a large range of parameter space and the computational 
cost involved in N-body codes proves to be a significant drawback. Furthermore, as was pointed out by Allgood et al. (2006), 
various groups differ in their methodology of determining halo shapes and hence differ in their results. Thus, the need for 
analytic investigations still remains. 

Triaxial collapse in an expanding universe has been under study for over four decades ( Icke| 1 973; White & Silk 1979; Nariai 
& Fujimoto 1972; Eisenstein & Locb 1995; Bond & Myers 1996. Lithwick & Dalai 20111. Various authors have considered 


it specifically in the context of axis ratio evolution. The seminal paper by Bardeen et al. (1986) used excursion set theory 
to compute the ellipticity and prolateness distribution of dark matter haloes forming under gravitational collapse of initial 
Gaussian random fields. Lee, Jing, & Suto (2005) combined the dynamics given by the Zeldovich approximation with the 
statistical properties of the linear field to obtain the probability distributions of axis ratios. |Sandvik et al.| ( |2007| ) and|pcsjacqucs 
120081 have used excursion sets to model the effect of the environment on halo ellipticity. 

Recent work by the authors (Nadkarni-Ghosh & Singhal [2014] hereafter paper I) analysed triaxial collapse in terms of 
the coupled dynamics of eigenvalues of three tensors: the Hessian of the gravitational potential, the tensor of derivatives of 
the velocity field and the deformation tensor. Paper I focussed on the joint dynamics of the gravitational and velocity fields. 
This paper focuses on the dynamics of the deformation tensor, which gives direct information about the evolution of the axis 
ratios. We start by reviewing the phase space equations and establish the notation (©• The details of the numerical runs 
are described next (©■ The shape of the triaxial object is quantified using three parameters: the minor to major axis ratio, 
s, the intermediate to major axis ratio, q, and the triaxiality parameter (Franx, Illingworth, & de Zeeuw 19911, T. These 
are defined in terms of the eigenvalues of the deformation tensor; the dynamics of the latter is completely determined by the 
phase space equations. S[4] gives the results and we discuss and conclude in S[5] 


2 THE PHASE SPACE EQUATIONS 


The physical system consists of a homogenous isolated ellipse evolving in a cosmological background consisting only of dark 
matter and dark energy with an equation of state w = — 1 (ACDM). In the absence of rotations, the evolution of the ellipse 
is completely determined by the three tensors: the Hessian of the gravitational potential, the tensor of velocity derivatives 
and the deformation tensor. Let Ad,i,A„,j and A a ,i {i=l,2,3}, denote the eigenvalues of these three tensors respectively. The 
dynamics of these nine variables is given by (see paper I for details) 
d\a,i 


— Au,i(l Aa,i) 


din a 

= — — [3f2 m (a)Ad,i — {f2 m (o) — 2flA(a) - 2} A V} i + 2A„ i J 


din a 

d\d,i 
din a 


»-(! + «) « + 


A d,i + g j A v,j 

' j= 1 


(la) 

(lb) 

(lc) 


+ ^A d,i + y ' (1 + A„,i) — ^ (1 + A V: i) 

{Adj — \d,i} ' {(1 — Ao,i) 2 (l + \v,i) — (1 - Aaj) 2 (l + A„,j)} 


+E 


where 


(1 - A a,i) 2 - (1 - Vi) 2 

* = !>.<• 


Let a be the scale factor of the background and m be the scale factors corresponding to the three axes of the ellipsoid. The 
three AiS and cqs are related as 


where 




OLi = 


roo 

dr 

L 
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&G 

Rf,Gauss (h Mpc) 

Mf , Gauss ( 10 14 M©) 

0.5 

14.26 

5.30 

0.6 

12.04 

3.18 

0.8 

9.05 

1.34 

0.9 

8 

0.91 

1 

7.12 

0.65 

1.2 

5.77 

0.34 


Table 1. Relation between the r.m.s. ctq and the mass scale Mf for a BBKS power spectrum with spectral index n s 

k 2 R 2 f 

2003) and a Gaussian filter W(kRf) = e 2 . The power spectrum was normalized such that ag = 0.9. 


1 < |Dodelson| 


and the evolution of Qs is given by 


^m(u) — 




TT 2 3 

<,,oVl 0 a 0 


H 2 a 3 


^a(o) = 1 — H m (a) 


(5) 


The subscripts ‘0’ denotes the value today. At early times, when the fluctuations are small, A d,i = A a ,i and A v ,i = —/(n m )A<;,i, 
where /(fl m ) is the linear growth rate, usually approximated as /(fi m ) « for a ACDM cosmology ( |Linder|[2005| ). The 
equations above were derived from the triaxial dynamics model of Bond & Myers| (|l996). 


3 NUMERICAL RUNS 


For Gaussian initial conditions, the initial distribution of the A^i is given by |Doroshkevich| |l970) 

15 3 ( 3 It 157 2 \ 


p( Ad,i, Ad, 2 , Ad, 3 ) = 


87 t \/5 o 


■ exp 


x (Ad,i — Ad, 2 )(Ad ,2 — Ad, 3 )(Ad,i — Ad, 3 ) 


( 6 ) 


’g ) 

where <jg is the r.m.s. fluctuation at the scale Rf, h = Ad,i + Ad ,2 + Ad ,3 and I 2 = Ad,iAd ,2 + Ad, 2 Ad ,3 + Ad,iAd, 3 . This 
distribution only assumes Gaussianity and does not depend on the details of the power spectrum. The power spectrum is 
relevant when relating og to a mass or radius scale: 

1 


°a{Rf,z) = 


(271 


J P[k, z)W 2 (kR f )d 3 k, 


(7) 


where P(k, z) is the linear power spectrum at redshift a and depends on the power spectrum at redshift zero through the linear 
growth factor (Dod elson] |20031 and W(kRf) is the window function. The mass scale Mf is related to smoothing scale Rf as 
Mf = An/3R 3 fpm, where pm is the homogenous background matter density. Table flj gives the relation between aa, radius and 
mass scales for the BBKS power spectrum (Bardeen et al. 1986[ ) with a Gaussian window function. We considered five different 
values of a g : 0.5,0.6,0.8,1 and 1.2, which covers a decade in mass scales. For each value of og five realizations were drawn, 
each consisting of 10 4 points. We focussed only on those initial points which correspond to overdense regions and are mildly 
non-linear today (a 0 = 1). This roughly corresponds to initial 8 values which lie within one standard deviation. At the initial 
time, chosen to be a ini t = 0.001, the system is linear and fi m « 1, implying A a ,i(ai n it) = \ v ,i( a init) = —A d,i{ai„it). 

The initial conditions were evolved from a = to a = 1 using the system of eqs. 0- The shape parameters of the ellipse 
at any intermediate time can be constructed from the A a s. 

The shape of a triaxial object is completely characterized by two parameters (a third parameter will correspond to setting 

the scale of the object). It is customary to use the ratios of smallest to largest and intermediate to largest axes defined as 

s = - - (8) 

u ’ma x 

Winter / n .\ 

q = - -, (9) 


where a m i n ,a ma;E and ai n t er are the smallest, largest and intermediate axis respectively. In terms of the A a parameters, these 
are 

1 A a.max 


S = 


Q = 


' a,min 


1- A a 

1 Aa inter 


1 A a,n 


The triaxiality parameter that characterizes the prolateness is defined as ( jFranx, Illingworth, fc de Zeeuw |l991 ) 


T = 


_ 2 -| _ 2 

^inter Q 


1 — S 


2 ' 


( 10 ) 

(11) 

( 12 ) 


Note that T is not an independent parameter, but derived from q and s. Often, in the literature, triaxiality is characterized 
in terms of the ellipticity e and prolateness p (Bardeen et al. 19861 and haloes are known to populate a triangular region in 
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fr=0.8, M=4.8x 10 14 M Q 


Scaled and shifted 




Figure 1. Probability density function of the minor to major axis ratio. The raw PDF is fit by a beta distribution and the scaled and 
shifted PDF represents has a universal fitting form given by the lognormal distribution. 


the e — p space (Porciani, Dekel, & Hoffman 2002; Dcsjacqucs 2008. Despali ,'Giocoli, fe Tormen||2014 1. These are equivalent 
variables, defined from the eigenvalues of the deformation tensor. While both sets are easy to measure in simulations, s and q 
have a more direct geometric interpretation. In this paper, we examine the evolution of s, q and T as a function of mass and 
redshift for a LCDM cosmology. Here we focus only on dark matter overdensities; the investigation of voids is left for future 


studies. We choose the cosmological parameters in accordance with the ACDM cosmology dictated by WMAP-7 (Komatsu 


et al.[2011): fi m,o — 0.29, fU,o — 0.71, Tts — 1} o '8 — 0.9, h — 0.73. The critical density of the universe reejuired to relate the 


Rf to mass values is p Ct o = 2.775 h 2 x 10 11 MgMpc 3 | Dodelson|2003 1. 


4 RESULTS 


Figure [I] shows the probability distribution function (PDF) of s. The left panel shows the variation with mass scale at a 
given redshift and middle panel shows the variation with redshift for a fixed mass. Smaller mass haloes (larger ctg) are more 
elliptical (smaller s) than the larger mass haloes. For a fixed mass scale, the halo ellipticity increases at lower redshifts. This is 


in agreement with the conclusion reached by Lee et al. (20051 which used the Zeldovich approximation to model the dynamics. 


Analytic work by |Bernardeau| ( |l994| | also suggests that rare peaks (more massive objects) are more spherical. However, these 
conclusions are contrary to the features observed in numerical simulations. This is as expected. Multiple numerical simulations 
(e.g. Jing fe Suto| [2002 Allgood et al.] |2006 1 have shown that accretion and merging play an important role in determining 
halo ellipticity. Low mass haloes collapse earlier and become more spherical due to accretion and mergers. Therefore, at a 
given redshift, the lower mass haloes are more spherical than higher mass haloes. Similarly, haloes of similar mass are more 
elliptical at a higher redshift since mergers are fewer at higher redshifts. The analytic model considered here deals with isolated 
ellipsoids; the effect of the environment, accretion and merging are not incorporated. Non-linear collapse simply amplifies the 
initial ellipticity and hence for any halo, ellipticity increases with decreasing redshift. 

Although analytical estimates and simulations are in contradiction, we find that a suitable transformation can reconcile 
the differences. In recent work, Bonamigo et al. 2014 (henceforth B14) analysed the Millennium XXL Simulations (Angulo 
et al. 2012) and showed that scaling the axis ratio s by an appropriate power of the peak height removed the dependence 


on mass and redshift resulting in a universal distribution. They considered cluster scale objects M ~ 10 M@ which roughly 
corresponds to the same range of ctg considered here (see Table [lj. Motivated by their results, we examined the effect of 
such a scaling on the parameter s. We find that pure scaling does not give rise to a universal profile; an additional offset is 
necessary. Define a new variable]]] 

s = si/ 0 ' 255 + s 0 (z,ctg), (13) 

where the peak height v = Here, 8i is the initial (linear) density of an individual point in the realization and ctg is the 
r.m.s. value. The offset parameter so(<jg,z) can be interpreted as a model for the physics missing in the analytic picture. 
Figure [2] shows the offset as a function of z for three different values of ctg. As expected, the offset is small at high redshifts 
and for larger mass scales (lower <tg). This function captures the fact that effect of mergers becomes increasingly important 
at late times. The scale at which mergers start to dominate also changes with time; lower mass scales are affected at earlier 
epochs because they collapse first. Fitting the offset as a function of ctg and z gives 

so(cr G , z) = A(ctg) log 10 ( App) + B(ctg), (14) 


l + z / 

where A(ctg) = 0.22 ctg + 0.44 and B(<jg) = —0.91 log 10 gg + 0.45. This fit is valid in the range ctg — 0.5 to gg — 1 and 


1 We differ in our definition of v from that of |Bonamigo et al.| < |2014| ) who use v = 5 c (z)/cr(M ), where 8 c (z) is the critical overdensity 
for collapse and a(M) is the r.m.s. fluctuation on the scale M. 
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Figure 2. Offset function so(z,a). Dots indicate the data points and the solid lines are the fit given by eq. (??). This function models 
the discrepancy between simulations and the analytic model presented here. 


the redshift range z = 1 to z = 0. The PDF of s is shown in the last panel of figure [I] It has a universal form given by the 
lognormal function. 

1 


p{s,n,a) = 


(In s — fi) 2 
6XP{ - 2tU } 


(15) 


Xy/2n a 

with fi = —0.49, a = 0.2, in agreement with that of B14. 

Figure [3] shows the PDF of intermediate to major axis ratio. The left and the middle panel show the PDF of the variable 
q. The definitions of q and s imply that q > s and also q < 1. Therefore, for the same mass scale and redshift, q has 
a narrower distribution than the corresponding s distribution in figure [l] The ratios s and q individually capture the two 
dimensional eccentricity along a cross-section perpendicular to the intermediate and minor axis respectively. To understand the 
full 3D geometry, it is generally the conditional probability p(q\s) which is the quantity of interest. Instead of the conditional 
probability, we consider the PDF of the s-dependent variable defined by [Schneider, Frenk, & Cole] ([2012]) 

< 16 > 

Their motivation for this transformation was that its range is the unit interval and hence it has the same support as the 
beta distribution. They considered the conditional probability p(q\s) and found s-dependent but redshift independent beta 
distributions. We do not consider the conditional probability, i.e., we do not bin in s space, instead we use the s value of the 
individual density peak (halo) and find that PDF of this transformed variable is independent of mass scale and redshift and 
also well fit by a beta distribution 

1 


p(q,a,P) = 


-X~\l-qY 


(17) 


B(a,0) 

with parameters a = 2.52 an (5 = 2.56 over the mass range a g = 0.5 — 1.2 and over the redshift range z = 99 to 2 = 0. The 
agreement at very early times, suggests that this scaling is encoded in the initial distribution, but the universality suggests 
that the non-linear evolution preserves it. To understand this universality better, we examined the pointwise evolution of 
q over the entire redshift range from the phase space dynamics. We focus on a single value of ctg and a single realization. 
For each point in the realization we examine the difference between the transformed ratio at any epoch q(a) and the ratio 
today qo = <?(a = 1)- Figure [4] shows the maximum difference over all points in the realization (for <jg = 0.5). This gives 
a conservative estimate of the difference for each epoch. We find that the difference stays small throughout the evolution. 
Thus, we can conclude that to within about a 1% error, q can be thought of as a invariant of the dynamics. The fact that the 
difference decreases monotonically is just a consequence of defining the differences with respect to the value today. Whether 
or not numerical simulations show a similar invariance/universality in the variable q remains to be checked. 

Figure [5] shows the distribution of the triaxiality parameter T. A perfectly prolate spheriod is one in which two axes are 
equal and less than the third (a m i n = ai n t e r < a ma x => q = s, T = 1) and a perfectly oblate spheriod is one where two axes 
are equal and greater than the third (a m i n < a,i n ter = Umax ==$• q = l,s < 1 ,T = 0). For general triaxial objects with no 
two axes equal, the range 0 < T < 1/3 is considered oblate, 2/3 < T < 1 is considered prolate. Haloes with 1/3 < T < 2/3 
are neither prolate nor oblate (Allgo od et al.|200 6 Schneider, Fren k, fe Cole|2012 |. The left panel of figure [H] shows that at a 
given redshift, lower mass haloes (higher og) are more prolate than higher masses and the right panel shows that prolateness 
is higher for lower redshifts. On the other hand, numerical simulations (Schneider, Frenk, fe Cole| |2012 | |Despali, Giocoli, fe| 
Tormcn 20141 and observations (Paz et ak||2011 ) suggest that most haloes are prolate, the more massive haloes being more 


prolate. Prolateness also increases with increasing redshift. These contradicting results again highlight the importance of post 
collapse processes, accounted for in simulations, in determining halo prolateness. 
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Scaled 



Figure 3. Probability density function of the intermediate to major axis ratio. Both the raw and the scaled PDFs follow the beta 
distribution. The latter is a universal profile with parameters a = 2.52 and 0 = 2.56. 



a 

Figure 4. q as an invariant of the dynamics. q(a) — qo gives a measure of the change in q for each point in the realization. The maximum 
over all points in the realization gives a conservative estimate of the overall change. The difference stays small throughout the evolution 
indicating that q may be treated as an invariant of the phase space dynamics. 


5 DISCUSSION AND CONCLUSION 


We have investigated the evolution of axis ratios using the dynamics of the eigenvalues of the deformation tensor. The 
eigenvalues are directly related to three parameters which quantify the shape of the triaxial object: the ratio of minor to 
major axis s, the ratio of intermediate to major axis q, and the triaxiality parameter T, which quantifies prolateness vs. 
oblateness. We find that, for a given redshift lower mass haloes are more elliptical and more prolate than higher mass haloes 
and the ellipticity and prolateness decreases with redshift. These trends are opposite to those predicted in simulations, but 
in agreement with other analytic studies such as|Lee et al. (2005]). This behavior is expected since it is known that mergers 
play a very crucial role in determining halo ellipticity. Furthermore, Zhao et al. (2003) have argued that more massive haloes 


undergo a faster accretion, which may account for partially compensating the trends predicted on the basis of isolated collapse, 
however, more recent work has suggested a mass independent accretion history (Ludlow et al. 2013). 

Inspite of this difference, we find that in transformed variables, we obtain a universal behavior for the two axial ratios. 
This is the main new result of this work. Motivated by the scaling proposed in B14, we find that the transformed variable 


cr=0.8, M=4.8 x10 14 M o 




Figure 5. Probability density function of the triaxiality parameter T. The tendency for prolateness increases with decreasing redshift 
and decreasing mass. 
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s = sv' 55 + so ( 2 , cr), where v is the peak height, has a universal log-normal profile. The main difference between the scaled 
variable s of B14 and ours is the additional offset so (z, a). This offset function can be thought to capture the missing physics of 
mergers and accretion and can be used to gain insights while building a more complete analytic model. The second transformed 
variable q = (q — s)/( 1 — s) proposed by Schneider, Frenk, & Cole (20121 was also found to have a universal PDF given by 
the beta distribution. By examining the point-wise variation of q, we find that the variation in q over the entire redshift 
range is very minimal indicating that q can be thought to be an invariant of the phase space dynamics. Similar phase space 
invariants have been studied in the context of density-velocity dynamics (Nadkarni-Ghosh 2013; Nadkarni-Ghosh & Singhal 


2014). Whether or not simulations suggest similar invariants remains to be checked. If such invariants exist, then they can 


be exploited to give additional constraints on cosmological parameters as was demonstrated based on spherical dynamics in 
Nadkarni-Ghosh]|2013j|. This is left for future investigation. 

The dynamical evolution equations used in this work were derived from the triaxial collapse model of of |Bond fc Myers] 
(1996). This model involves many approximations. The ellipsoid is rotationless and isolated and the effect of external tidal 
forces is modeled in terms of the axes of the ellipsoid. In order to make more realistic predictions, it will be necessary to 
consider models such as the one proposed by Desjacques (2008) which more accurately takes into account the effect of the 
environment and/or employ the formalism developed by Nariai fe Fujimoto| ( |1972| and Eiscnstci n fe Loeb| fl995| ) which allows 
for rotations. The main motivation to build increasingly accurate analytic models, inspite of the increasing accuracy of N-body 
simulations, comes from the point of view of precision cosmology. Given the range of unconstrained cosmological parameters 
and the increasing size of data sets, computing time becomes a limiting factor. Analytic alternatives are hence necessary and 
studies like the one presented here pave the way for more sophisticated model building in the future. 
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